%%傅里叶变换加水印源代码
%% 运行环境Matlab2021a

clc;clear;close all;
alpha = 20;%水印的权重 可以调节
%% 读取原图和水印图片

im = double(imread('football.jpg'))/255;
mark = double(imread('watermark.png'))/255;
figure, imshow(im),title('Original image');
figure, imshow(mark),title('watermark');
bw1=im2bw(mark);
figure, imshow(bw1);
%% random

imsize = size(im);
TH=zeros(imsize(1)*0.5,imsize(2),imsize(3));
TH1 = TH;
TH1(1:size(mark,1),1:size(mark,2),:) = mark;
M1=randperm(0.5*imsize(1));
M=sort(M1);
N1=randperm(imsize(2));
N=sort(N1);
save('encode.mat','M','N');
for i=1:imsize(1)*0.5
    for j=1:imsize(2)
        TH(i,j,:)=TH1(M(i),N(j),:);
    end
end
mark_ = zeros(imsize(1),imsize(2),imsize(3));
mark_(1:imsize(1)*0.5,1:imsize(2),:)=TH;
for i=1:imsize(1)*0.5
    for j=1:imsize(2)
        mark_(imsize(1)+1-i,imsize(2)+1-j,:)=TH(i,j,:);
    end
end
figure,imshow(mark_),title('encoded watermark');
%% Add watermark

FA=fft2(im);%给原图做傅里叶变换
figure,imshow(FA);title('Spectrum of original image');
FB=FA+alpha*double(mark_);%给图像加上水印
figure,imshow(FB); title('Spectrum of watermarked image');
FAO=ifft2(FB);
figure,imshow(FAO); title('Watermarked image');
%% 涂抹攻击

attack = 1-double(imread('attack.jpg'))/255;
figure,imshow(attack);title('涂抹攻击');
FAO_ = FAO;
attsize = size(attack);
for i=1:attsize(1)
    for j=1:attsize(2)
        if attack(i,j,1)+attack(i,j,2)+attack(i,j,3)>0.3
            FAO_(i,j,:) = attack(i,j,:);
        end
    end
end
figure,imshow(FAO_);
%解出水印
FA2=fft2(FAO_);
G=(FA2-FA)*2;
GG=G;
figure,imshow(G);title('extracted watermark');
for i=1:imsize(1)*0.5
    for j=1:imsize(2)
        GG(imsize(1)+1-i,imsize(2)+1-j,:)=GG(i,j,:);
    end
end
figure,imshow(GG);title('extracted watermark');
%% attack by cutting

s2 = 0.8;
FAO_ = FAO;
FAO_(:,s2*imsize(2)+1:imsize(2),:) = FAO_(:,1:int32((1-s2)*imsize(2)),:);
figure,imshow(FAO_);title('剪贴后的图像');
%解出水印
FA2=fft2(FAO_);
G=(FA2-FA)*2;
GG=G;
for i=1:imsize(1)*0.5
    for j=1:imsize(2)
        GG(imsize(1)+1-i,imsize(2)+1-j,:)=GG(i,j,:);
    end
end
figure,imshow(GG);title('extracted watermark');
%% 旋转

FAO_ = FAO;
I1=myimrotate(FAO_,90);     %调用myimrotate()函数旋转30° 
I2=myimrotate(im,90);
figure,imshow(I1);title('旋转后的水印图像');
figure,imshow(I2);title('旋转后的原图像');
FA=fft2(I2);
FA2=fft2(I1);
figure,imshow(FA);title('原图旋转后的傅里叶谱');
figure,imshow(FA2);title('水印图旋转后的傅里叶谱');
J = imrotate(FAO_,90);
figure,imshow(J);title('旋转后的水印图像');
FJ=fft2(J);
figure,imshow(FJ);title('原图旋转后的傅里叶谱');
J1 = imrotate(im,90);
figure,imshow(J1);title('旋转后的水印图像');
FJ1=fft2(J1);
figure,imshow(FJ1);title('原图旋转后的傅里叶谱');

G=(FJ-FJ1)*2;
abs(fftshift(G));
figure,imshow(G);

G=-(FA2-FA)*2;
figure,imshow(G);
%% 椒盐噪声

FAO_ = imnoise(FAO,'salt & pepper',0.01);
figure,imshow(FAO_);
FA=fft2(im);
FA2=fft2(FAO_);
G=(FA2-FA)*2;
GG=G;
for i=1:imsize(1)*0.5
    for j=1:imsize(2)
        GG(imsize(1)+1-i,imsize(2)+1-j,:)=GG(i,j,:);
    end
end
figure,imshow(GG);title('extracted watermark');
%%
function [ A ] = myimrotate(B,degree)                                 %定义旋转函数，degree为要旋转的角度
[r,c,d]=size(B);                                                      %获取输入图像B的行r、列c和通道数d,为了旋转彩色图像所以有必要得到通道数d
nH=round(r*abs(cosd(degree))+c*abs(sind(degree)));                    %旋转图像后得到的新高度，“round()函数四舍五入“
nW=round(c*abs(cosd(degree))+r*abs(sind(degree)));                    %旋转图像后得到的新宽度
A=zeros(nH,nW,d);                                                     %定义生成目标图像的行列以及通道数
M1=[1 0 0;0 1 0;-0.5*nW -0.5*nH 1 ];                                  %坐标系变换矩阵M1
M2=[cosd(degree) -sind(degree) 0;sind(degree) cosd(degree) 0;0 0 1];  %角度旋转变换矩阵M2，顺时针方向
M3=[1 0 0;0 1 0;0.5*c 0.5*r 1];                                       %坐标系变换矩阵M3
    for i=1:nW
        for j=1:nH
            temp=[i j 1]*M1*M2*M3;                                    %得到旋转后的矩阵temp
            y=temp(1,2);                                              %y取矩阵temp的第一行第二列,y对应j，为高度
            x=temp(1,1);                                              %x取矩阵temp的第一行第一列,x对应i，为宽度
            y=round(y);                                               %y四舍五入取整
            x=round(x);                                               %x四舍五入取整
           if(x>=1&&x<=c)&&(y>=1&&y<=r)                               %判断的得到的(x,y)点是否在原图像上
               A(j,i,:)=B(y,x,:);                                     %将原图像的像素点赋值给对应的旋转后图像上的点
           end                                                        %i,x对应的是列，即宽，而j,y对应的是行，即高，以x为横坐标，y为竖向纵坐标
        end
    end
end